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A mean-field Fokker-Planck equation approach to the dynamics of ferrofluids in the presence of a magnetic 
field and velocity gradients is proposed that incorporates magnetic dipole-dipole interactions of the colloidal 
particles. The model allows to study the combined effect of a magnetic field and dipolar interactions on the vis- 
cosity of the ferrofluid. It is found that dipolar interactions lead to additional non-Newtonian contributions to the 
stress tensor, which modify the behavior of the non-interacting system. The predictions of the present model are 
in qualitative agreement with experimental results, such as presence of normal stress differences, enhancement 
and different anisotropy of magnetoviscous effect and the dependence of the viscosity on the hydrodynamic 
volume fraction. A quantitative comparison of the concentration dependence of the magnetoviscosity shows 
good agreement with experimental results for low concentrations. 

PACS numbers: 47.65.+a Magnetohydrodynamics and electrohydrodynamics; 75. 50. Mm Magnetic liquids; 47.32.-y Rota- 
tional flow and vorticity; 83.80.Hj Suspensions, dispersions, pastes, slurries, colloids; 05.20.Dd Kinetic theory 



I. INTRODUCTION 

The so-called magnetoviscous effect in ferrofluids - the dependence of the rheological behavior of colloidal suspensions of 
nano-sized ferromagnetic particles in a carrier liquid on external magnetic fields - has attracted quite some attention in the 
recent literature, both from a scientific and an application point of view GIB 0] ■ In dilute ferrofluids, the magnetoviscous 
effect is successfully explained by the hindrance of rotation of individual, non-interacting magnetic dipoles Q]. Experiments 
on commercial ferrofluids have revealed quantitative and qualitative different behavior compared to the dilute regime, such as 
an order of magnitude increase of the magnetoviscous effect, a dependence not only on the local vorticity but also on the local 
strain rate of the flow, shear thinning behavior and the occurrence of normal stresses The failure of the kinetic model 

0] to describe these phenomena is generally attributed to the neglect of magnetic dipolar interactions of the colloidal particles. 
Since dipolar interactions are long-ranged, they become important at concentrations as low as a few percent which are easily 
attained in commercial ferrofluids. The present contribution provides an extension of the non-interacting model to moderately 
concentrated ferrofluids with weak dipolar interactions. 

A general statistical theory of magnetic fluids that covers dilute as well as concentrated suspensions with arbitrary strength of 
dipolar interactions is not available at present. For strong dipolar interactions, A > 1, the formation of chain-like aggregates is 
expected. A corresponding kinetic model has been proposed in Ref. |(J. In many ferrofluids, the interaction parameter A is of 
order unity at room temperature, implying that thermal energy is sufficiently strong to prevent permanent aggregation |3|. For 
example, rheological measurements on magnetite based ferrofluids reported in Ref. [7fl estimated A to be 0.2. We here propose 
an extension of the kinetic model 0|, that does not assume the existence of permanent chainlike aggregates but that incorporates 
magnetic dipole-dipole and excluded volume interactions in a mean-field approximation. The model is expected to apply in the 
dilute and semi-dilute, weakly interacting regime, A < 1. Several results are obtained: A correction to the Langevin function 
describing the equilibrium magnetization of non-interacting magnetic dipoles is obtained in agreement with results of Ref. Q. 
The additional contributions of dipolar interactions to the hydrodynamic stress tensor are worked out. Contrary to the case of 
non-interacting magnetic dipoles, the stress tensor now depends also on the symmetric part of the velocity gradient. In case of 
uniaxial symmetry, the hydrodynamic stress tensor is of the same form as in the Ericksen-Leslie theory of nematic liquid crystals. 
The predictions of the present model are in qualitative agreement with experimental results and address the shortcomings of the 
non-interacting model 0|, such as the modified concentration dependence and anisotropy of the magnetoviscosity and the 
presence of normal stresses. Some quantitative comparison to the experimental results of Ref. 0] are offered also. It should be 
mentioned, that a similar approach has been proposed already in Ref. |9||. In |9[, however, no dependence on the rate-of-strain 
tensor is considered, while such a dependence has been found experimentally in Ref. llfJl . In addition, we keep higher order 
contributions in A compared to Ref. |9|] ■ 

This paper is organized as follows: In Sec. [H] the kinetic model of semi-dilute ferrofluids is developed for equilibrium 
conditions. It is shown that the equilibrium behavior of this model is agrees with the results of Ref. J8|. In Sec. the 
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kinetic model is extended to describe the dynamics of ferrofiuids in the presence of an external flow field. The kinetic model 
is supplemented by the definition of the hydrodynamic stress tensor. The rheological behavior of the present model is studied 
in Sec. IIVI In case of uniaxial symmetry, the hydrodynamic stress tensor is found to be of the general form proposed in the 
Ericksen-Leslie theory of nematic liquid crystals. Explicity expressions for the viscosity coefficients are obtained in case of 
weak flow and compared to experimental results. Finally, some conclusions are drawn in SecfVll 



n. MODEL DEFINITION AND EQUILIBRIUM PROPERTIES 

Consider a system of N interacting, spherical colloidal particles in a volume V. All particles are assumed to be identical, 
ferromagnetic monodomain particles of diameter d. Let Xi = {r'*', uW} denote the five-dimensional vector describing the 
position j»W and orientation u'*' of particle i, u^ 1 ' ■ u'*' = 1. The particles are assumed to carry a permanent magnetic moment 
= iavM\ The total interaction potential may be written as 

N 

2—1 i<j i<j 

The first term denotes the potential energy of an ideal paramagnetic gas in the local magnetic field H\ oc . The second term is the 
potential energy of the non-magnetic system, where = w s (rij), with = r« - r®, t 2 j = rf,-, is a spherical symmetric, 
short range, repulsive potential. In particular, we consider the case of hard spheres, w s (r) = oo if r < d and zero otherwise, and 
soft spheres, w s (r) = (r/d)~ 12 . The energy of two magnetic dipoles is described by 



f3wt?{r l0 ,u^M 3) ) = -3A(d/r 4J ) 3 « (l) ■ r^a ■ (2) 

where = r^r,^ and a" = (a + a T )/2 — (tr a)l/3 denotes the symmetric traceless part of the matrix a. The dimensionless 
interaction parameter 

A = AirLi k B Td 3 (3) 

is given by the ratio of the magnetic dipole-dipole energy of two colloidal particles of diameter d in contact over the thermal 
energy. It is well-known that due to the long range nature of the dipolar interactions, the magnetic properties of the system 
depend on the geometry of the container. In order to deal with this situation, we follow Ref. |$j and introduce a virtual cut-off 
radius R c of the dipolar interactions. The effect of dipole j on dipole i with r,j > R c is treated in a continuum approximation. 
Within the Weiss model, the resulting far-field contribution leads to a local magnetic field H\ oc which is given by 

ifioc = H + ^M, (4) 

where M denotes the magnetization of the sample. Finally, the virtual cut-off is removed, R c — > oo. For further details see 
Ref. 

Exact results for the thermodynamic properties of the model system Q are not available. Since the typical volume fraction of 
magnetic material in ferrofiuids is low, the free energy of the system is conveniently expressed by the virial expansion. Let f(u) 
denote the one-particle distribution function of finding the orientation u of an individual colloidal particle. The normalization 
is chosen such that Jd 2 uf(u) = 1, where integration over the three-dimensional unit sphere is denoted by Jd 2 u. The free 
energy functional per particle of the system, F [/] , may be split into an ideal, Fq [/] , and an excess part, F ox [/] . The ideal system 
consists of an ideal gas with number density n = N/V and an ideal paramagnetic gas, 



F [f] = k B T 



Inn — 1 — J d 2 uf(u)u ■ h\ oc + / d 2 uf(u) In f(u) 



(5) 



where h\ oc = [iH\ oc / k B T denotes the dimensionless local magnetic field. Boltzmann's constant and temperature are denoted 
by k B and T, respectively. For low concentrations, the excess part may be written in terms of the virial expansion as fllll 



FeM] = ~^nk B T Jd 2 ujd 2 uif(u)f(ui)b 2 (u,ui) + 0(n 2 ). (6) 

The function 62 is defined by 

b 2 (u,u') = Jd 3 r(e-^ r ^-l) 9sp (r), (7) 
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where g sp denotes the pair correlation function of the reference system (wf^ — 0) and = k^T. The function 62 can 
be interpreted as the change of the second virial coefficient due to the dipolar interactions. The integration over the three- 
dimensional spherical volume is denoted by Jd 3 r. 

The function b 2 , Eq. 0, can be expressed as a power series in the interaction parameter A. Using the fact that the pair 
correlation function g sp {r) is spherical symmetric, we arrive at 

6 2 ( u X) = 24«]TA fc -— C2jfc .„')* (8) 

fc=2 ^ ' " J n 

where v = nd 3 /6 denotes the hydrodynamic volume of the colloidal particles. The numerical coefficients 

/•OO 

C2,* = 3(fc-1) / dxx 2 - 3k g(x). (9) 
Jo 

depend on the particular choice of the short range potential w s and corresponding pair correlation function. Due to the spherical 
integration volume the term k = 1 is missing in the sum of Eq. (|8|. In the low density limit, g sp {x) can be approximated 
by g sp (x) « cxp(— f3w s ). If we consider the case of hard spheres, g sp (x) can thus be identified in the low density limit with 
the Heaviside step function at x = r/d = 1, In this case we have c 2 ,k = 1. For power law repulsions, (3w s (x) — x~ v , the 
coefficients c 2tk are given by c 2 , k = fcr(fc), where k = 3(fc — \)/v. The so-called soft sphere potential is recovered for v = 12. 
In this case, the coefficients C2,fc are close to 1 for k < 5 and therefore similar to the value of the hard sphere system. 
Inserting the expansion (|SJl into Eq. © one obtains 

00 

F[f] = Fo[f] - q>k B TY,* k C2,kG k [f} + 0(<j> 2 ). (10) 

fc=2 

We have found explicit expressions of functionals G k for k < 5, 

2 / . , , , 10\ 

G 2 [/] = - (uu>:(«u) + — (11) 
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G 3 [f] = -— \(u a upv*y ){u a u u 1 ) - • ( u ) J (12) 



Gi[f] = ({UaUpU^US ){u a U[ 3 U~ 1 U S ) + y( tt «): <uu> + y ) . (13) 



Angular averages of arbitrary functions a(ti) with respect to the distribution function / are denoted by 

(a) = / d 2 u a(u)f(u). (14) 



Note, that functionals G k [/] depend on the distribution function only via moments of / up to order k. 

The equilibrium distribution / cq is found by minimizing the functional dlOt subject to the constraint of fixed normalization, 
/ d 2 u / C q(u) = 1. The result reads f C q(u) = cxp [— (3V MF (u; f eq )]/Z eq , where Z eq denotes the normalization constant. The 
mean-field potential is 

f3V MF (u; /) = -u • h loc d>f] A fc c 2lfc ^i. (15) 

^2 

Note that due to the occurrence of moments in Eq. Jl 5i the equilibrium distribution / eq has to be determined self-consistently 
signaling the mean-field character of the present model. Linearization in the volume fraction <fi leads to 



/cqO) = /«.(«) 



1 + — \ 2 <\>c 2 , 2 L 2 (a s ) [P 2 {u ■ H) - L 2 (a s )} 
A 3 0c 2 ,3 (L 3 (a s ) [P 3 (u ■ H) - L 3 (a s )} - 21L 1 (a a ) [u -H- Li(a s )}) 
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(16) 



\Uc2A (L 4 (a a ) [Pi(u ■ H) - L 4 (a s )} + 20L 2 (a s ) {p 2 (m • H) - L 2 (a s )} 



3675 
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where dots denote higher order terms in A. In Eq. d!6i . we have used functions Lj which are defined recursively by Lj+i (x) = 
Lj_i(x) - (2j + l)Lj(x)/x with L (x) = 1. 
The functions 

^ = A ^ Ujl (17) 

47r smn(a) 

are the equilibrium distribution functions in the absence of dipolar interactions. In Eq. illl we have introduced the Langevin 
parameter a, h = fiH/k-gT = aH with H the unit vector in the direction of the magnetic field. The macroscopic magnetization 
is expressed as M = M sat (u). Thus, the local field h \ oc can be expressed as h.i oc = a s H, with the effective Langevin parameter 
a s = a + x lj Li(a s ). The Langevin function is defined by L\{x) = coth(x) — x^ 1 and the Langevin susceptibility is x L = 8<^A. 

Evaluating the free energy functional fllUI > with the equilibrium distribution / oq one obtains the equilibrium free energy 
F(a) = F[/ eq ] up to linear order in <f>, 

F(a s )/k B T = In ( ^ j - ^ A fe c 2 , fc G fe (a s ). (18) 



k=2 



Functions Gk{a s ) are defined by Gk{a s ) = Gk[f cq }- Explicit expressions for the first functions Gk, obtained from Eqs. (II 11 - 
I13> combined with d!6> . are given in the appendix|X] We have confirmed that Eqs. dl 81 and JA1HA3> agree with the results of 
Ref. J3l in the case of hard spheres where c 2 ,k = L The advantage of the present formulation compared to the results of Ref. @] 
is, that Eqs. (IA1HA3> simplify the discussion of the properties and asymptotic behavior of the functions Gi(x). 

Define equilibrium order parameters by 5* q = (Pj (u ■ H)) cq , where Pj are Legendre polynomials of degree j. The function 
Sj q can be obtained by multiplying Eq. dl 61 by u ■ H and subsequent integration over u, or from Eq. dl 81 by S^ q (a s ) = 
dF(a s ) j da s . Linearization in the small quantity <fi leads to the final result 

00 

S^ia) = L^a) + xMa)^) + cf>Y / C2,k\ k G> k (a), (19) 

fc=2 

where prime denotes the total derivative. Eq. dl 91 is identical with Eq. (4.24a) of Ref. HI for C2,fc = 1. Fig. [2 shows a 
comparison of Eq. ( I19> with C2,k = 1 for k < 4 and C2 : k = for k > 4 to the results of molecular dynamics simulations 
of Ref. 1H2I1 . The volume fraction was chosen as (f> = 0.157 and the dipolar interaction was A = 1, which is rather large for 
the present study. From Fig.[^one notices that appreciable corrections to the Langevin magnetization occur for intermediate 
values of a. For the parameters considered, the simulation results are in remarkable agreement with Eq. d!9t if truncated at 
k = 4. In the regime of strong dipolar couplings, A > 2, where corrections to the Langevin magnetization are stronger fl2il . 
the truncation of the expansion d!9t is not admissible. This regime is left for future research. The truncation of the expansion 
at k = 2 gives similar results to the truncation at k = 4 for A ~ 1. In Fig.^ we include also the function Li(a s ), which is 
obtained by neglecting all higher order corrections in A in Eq. d!9t if the linearization in <f> is not performed. As has been noted 
in Ref. 11211 . this approximation describes their numerical data very well for the present choice of parameters. For later use, we 
give the expression for in linear order in <j>, 

00 

S 2 °» - L 2 (a) + Xl L 2 (c»0Xi(a) + 0^ c 2 . k X k J' k (a), (20) 

fe=2 

where 

= -^L 2 {a)[18L 4 (a) + 10L 2 (a) + 7 - 35L 2 (a) 2 ]. (21) 
525 



III. MEAN-FIELD KINETIC MODEL 

The model introduced in Sec.[n]is now extended to describe the nonequilibrium dynamics of ferrofluids in the presence of 
an external flow field v(r). The one-particle distribution function /(it) now becomes time-dependent, f(u;t), and denotes 
the probability density of finding the orientation u of an individual colloidal particle at time t. For convenience, the explicit 
dependence of / on t is frequently suppressed in the sequel. The orientational dynamics of a ferromagnetic colloidal particle 
under the combined action of the local vorticity of the flow fl = \ V x v, Brownian motion, and the action of the potential V 
is given by the kinetic equation 1301 

d t f - C ■ [Clf] + C ■ D Iot [Cf + fC((3V cS )] . (22) 
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The rotational diffusion coefficient is D Iot = l/(2r rot ), r lot = Srjsv/knT denotes the rotational relaxation time. The rotational 
operator is C = u x d/du with d/du the gradient on the unit sphere. In the absence of flow, we assume that the effective 
potential V cfi can be identified with the static mean-field potential V MF , Eq. H5\ . A similar approach was proposed by one of 
the authors in Ref. fl3ll in order to describe the dynamics of nematic liquid crystals. For simple fluids such an approach has been 
proposed and tested recently in Ref. Q. 

In the presence of a symmetric velocity gradient D = Wv , an additional contribution to the effective potential V cS of 
the kinetic equation d22b arises. In the case of non-spherical particles, this contribution leads to the so-called flow alignment 
phenomenon I13tfl5ll . In the present case, the additional contribution is due to flow-induced structures that can be formed even 
in a hard sphere system. The distortion of the pair correlation function due to shear flow has been studied experimentally fl^l 
and theoretically IPl . For small distortions, the pair correlation function g(r; t) satisfies the time evolution equation fl7ll 

d t g + r ■ (V r u) ■ V r .g + -(g - g sp ) = 0, (23) 

T 



where r denotes a translational relaxation time. The stationary solution to Eq. d23i is given by 

g(r) = g sp (r) - tD: ffrg' sp {r). (24) 

Results of recent nonequilibrium molecular dynamics simulations of a simple dipolar fluid confirm that Eq. J24I provides a 
reasonable description of the shear-induced distortion of the pair correlation function lisll . The distortion of the pair correlation 
function leads to an additional contribution to the effective potential which to lowest order in A reads 

V D (u;f) =n I d 2 u' /(«') / d 3 rwfi(r,u,u'){g(r) - g sp (r)). (25) 



Eq. J25l > has an immediate interpretation as the flow-induced modification of the (self-consistently averaged) mean dipolar 
interaction potential. The effective potential V ctl entering the kinetic equation J22b is obtained as V ctl = V MF + V u . Inserting 
Eq. &241 into d25i . the kinetic equation J22I takes to form 

dtf = C ■ [{fl - a oX ^u xD-(u)- D r C(pV MF )}f] + D r C 2 f, (26) 

where (Tq = 3r/(5r ro t)- In Eq. J26L we have assumed g sp {r) — > 1 for r — > oo and g sp (0) = due to excluded volume 
interactions. The kinetic equation for non-interacting dipoles, Ref. @], is recovered from Eqs. i26\ in the limit A — > 0. The use 
of Eq. d24b for the flow-induced distortion of the pair correlation function limits the validity of Eq. i26\ to weak flows. More 
precisely, we expect Eq. d24b to apply for t\D\ < 1. Since in Eq. d25t we have kept only the lowest order term in A, further 
considerations are restricted to the regime of weak dipolar interactions A <C 1. In this regime, the mean-field potential dl 5I > 
simplifies to 

f3V MF (u:f) = -u ■ h loc - ^X 2 cj)C2,2^T : (mi). (27) 

5 

As has been noted in Sec. |n] the truncation 1271 of the mean-field potential provides a good approximation for the equilibrium 
magnetization even for x L w 1, A « 1, 

It is interesting to compare the hydrodynamic drag appearing in Eq. d22i to the corresponding terms in case of ellipsoidal- 
shaped particles with axis ratio r. For the latter, the contribution of D to the hydrodynamic drag is quadratic in u and proportional 
to the so-called shape factor B = (r 2 — l)/(r 2 + 1) as well as to the rotational relaxation time r rot lfl5ll . Spherical particles 
correspond to r = 1 and B = 0. If it would be allowed to replace u(u) in Eq. d26t by uu, the flow contribution of the present 
model would be identical to a dilute solution of ferromagnetic ellipsoidal particles with shape factor B = —<tqx l - A negative 
value of B is characteristic of an oblate ellipsoid. This finding has an intuitive interpretation since dipolar interactions favor a 
head-to-tail arrangement in contrast to side-side configurations. Thus, the particles effectively appear extended in the direction 
perpendicular to the magnetic moment. In Sec. llVI the validity and limitations of the analogy are discussed for special quantities. 

For irrotational flows Q = 0, the stationary solution to the kinetic equation i22\ reads / st = exp (— (3V e )/Z st , where Z st 
denotes the normalization constant. Note, that / st is formally identical to the equilibrium distribution / eq , if the magnetic field 
hloc is replaced by the effective field h c g = h\ oc — (6/5)x l t£) • (u). 

From the kinetic equation i22\ a hierarchy of moment equations can be derived. Multiplying Eq. i22\ by u and integration 
over u leads to 

d t (u) = Q x (u) + oox L i( uu ) ■ D ■ (u) — D ■ («)) + D mt {h - (uu) ■ h) 

+ D rot x L {(u) - (uu) ■ (u)) + ic 2 ,2-DrotX L A((mr) • (u) - (uuu): (uu)) - 2D Tat {u). (28) 

5 
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Analogously, higher order moment equations are derived. In order to study the magnetoviscous effect, the kinetic equation d22> 
has to be supplemented by the expression for the hydrodynamic stress tensor T. Considering only rotational degrees of freedom, 
the stress tensor is antisymmetric and has the form T lot = 3r) s (f>e ■ (tt — (w p )), where e is the conventional total antisymmetric 
(Levi-Civita) tensor of rank three. The average angular velocity (uj p ) of the colloidal particles can be inferred from the kinetic 
equation (I22> which can be cast into the form dtf = C ■ [<*> p /]. Thus, T rot is given by 2T lot = ne ■ {CV cS ) and becomes 
upon inserting Eqs. (125 \ and (I27> 



-,rot nk BT ( _ ri .-V/..\ / . _ \ T -n 6 



T rot = ^ a[H{u) _ {u)H] _ TxAD . {u){u) {u){u) , D] y (29) 

One has x L = in the non-interacting model fl, so that the hydrodynamic stress arises from hindered rotations of individual 
magnetic moments in the magnetic field. The configurational contribution to the hydrodynamic stress is given by 



^jconf 



y J d 3 rj du J d«7(«)/(t»')fl (a) (r J t»,« / )rVr*ia(r > «,«'), (30) 



where $12(1", u, u') = w sp (r) + wf^if, u, u') is the full two particle interaction potential and g^ denotes the full pair cor- 
relation function. The distortion of from its equilibrium giq 1 is assumed to be described also by Eq. (I23> where g is now 
replaced by g^ 2 \ In this case, the stationary solution reads to leading order 

g^ ( r> u , u') = gW(r, u, u') - r(V r v) : rV r ^» (r, u, u 1 ). (31) 

Note, that Eq. (13 1 1 reduces to Eq. 1241 only if the (anisotropic) contribution of the dipolar potential to the pair correlation function 
is neglected. 

Inserting Eq. J3 1 1 into (I30i . the deviation of T con from the equilibrium stress becomes 



with the viscosity tensor 



conf 



2k B T 



a rriconf conf x~? /oon 
AT ^ = %u a p"aV^ (32) 



J d 3 r JduJ du' /(ti)/(« / )ffg ) (r,« J ti')r p (V 1/ #i2)r a (V i 8*i 2 ). (33) 



From Eq. J33i we notice that the viscosity tensor j^ conf obeys rf°j£p = ff^L, and is positive semi-definite, i.e. V c l ^a/3 a ^ a ai3 ^ 
for arbitrary second-rank tensors a. 

Inserting the spherical and the dipole-dipole interaction potential into Eq. (I33i yields symmetric and antisymmetric contribu- 
tions to the stress tensor, such that the total hydrodynamic stress tensor becomes 



with the symmetric part, 



2m - | (ci - 3)a{u)A D - 7 -a [W ■ (u) (u) - (u) (u> ■W]+a(c 1 - 3) [D ■ (u) (u) + (u> (u> • D] (35) 



and the antisymmetric part 



^(H(u)-(u)H). (36) 



In Eq. J35I we have introduced W = [(\7v) T — Vv]/2 and the quantity a = 2r] s aX<p 2 , where a = 72r/(35r ro t). The shear 
viscosity of the isotropic suspension is 

=%(l + ^ + 60 2 ), (37) 

where i] s is the shear viscosity of the pure solvent, b = \c^a. Note, that the isotropic interactions alter the value of the Newtonian 
viscosity while dipolar interactions lead to additional, non-Newtonian contributions to the stress tensor. Diagonal contributions 
to T s have been suppressed in Eq. (1351 since they can be compensated by the scalar pressure. Note also, that the stress tensor T 
is symmetric in the absence of an applied magnetic field. Quantities c\ and C4 result from the short range interaction contribution 
to the stress tensor and depend on the detailed form of the interaction potential, 

/>oo 

c k = / dxx k [l3w' a (x)} 2 g sp (x), (38) 
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where prime denotes the total derivative and x = r/d. The coefficients Cfc are non-negative, Ck > 0. Evaluation of the 
coefficients c\ and C4 for the case of hard spheres suffers from the discontinuity of the potential. For power law repulsion, 
f3w s (x) — x~ h ', the integration in J38I can be done analytically in the low density regime, giving = vT((2v + 1 — k)/v). 
For soft spheres, v = 12, we find c\ = 12 and C4 w 11.0. We mention that the tp 2 contribution to r]o in Eq. I37i stems from 
the non-magnetic interactions of the colloidal particles. This contribution is of the form G s r, where G s is the Born-Green 
expression of the equilibrium shear modulus of a system of spherical particles interacting with the potential w s 0. 

Eqs. 1261 1. (27\ , 1341 . ( I35l l and d36l constitute the mean-field kinetic model of the nonequilibrium dynamics of dilute, weakly 
interacting ferrofiuids proposed here. The present model extends the work of Ref. |4|] to the regime A *C 1. If the dimensionless 
dipolar interaction parameter A goes to zero, the model of non-interacting magnetic dipoles proposed in Ref. 0] is recovered 
with a renormalized zero-field viscosity t]q. 

The present model contains the additional parameter a which is a measure for the ratio of translational and rotational relaxation 
times. Simple estimates of the translational and rotational relaxation times lead to a w 6(ro/d) 2 , where ro is a typical length 
scale associated with the formation of flow-induced structures. If ro is identified with the typical distance between colloidal 
particles, the parameter a can be estimated as a ps 20~ 2 / 3 . On the other hand, if ro is identified with the diameter d of the 
colloidal particles one obtains a = 6. These estimates of a agree with each other for <j> « 0.2, which is slightly above the 
expected range of validity of the present model. Alternatively, if one requires the expression d37i for the zero-field viscosity 
to correspond to Batchelor's result, Eq. (I37> with 6 = 6.2 (see Ref. 0] and discussion therein), the parameter a is given by 
(7 = 186/(35c4) and thus related to the interaction potential w s . In case of soft spheres, agreement with Batchelor's result leads 
to a « 0.5. In the sequel, we consider a as parameter with a = 0(1). Fig.[2]shows zero-field zero-shear viscosity 770 of a 
kerosene based ferrofluid as a function of the volume fraction </>. The data are taken from Ref. 0. From Fig.|2]we notice that 
Batchelor's result describes the experimental data well for volume fractions <f> < 0.25. 



IV. RESULTS FOR UNIAXIAL SYMMETRY 

The stress tensor T, Eq. (I35I36> . depends explicitly on the first moment of the distribution function only. However, all the 
moments are coupled dynamically as can be seen already from the first moment equation d28l . such that the values of second and 
third moments are needed in order to determine the first moment and the stress tensor. Therefore, closed form equations for the 
dynamics of the stress tensor d34 | in terms of low order moments necessarily introduce approximations to the underlying kinetic 
model. In a previous work 12011 . we have studied the assumption of uniaxial symmetry of the distribution function applied to 
the non-interacting kinetic model of Ref. 01. Motivated by the good accuracy of the assumption of uniaxial symmetry for that 
model found in Ref. we employ this assumption also for the present case. 

In the uniaxial phase, the distribution function f(u; t) is symmetric with respect to rotations around the director n, f(u; t) = 
funi(u ■ n; t), such that / can be represented as 

f uni (u ■ n; t) = i- £ —L-SjMPjiu ■ n). (39) 

3=0 3 

The scalar orientational order parameters Sj are defined as Sj = (Pj{u ■ n)), where, as before, Pj are Legendre polynomials. 
In case of uniaxial symmetry, the first moments take the form 

(u) = Sin, (uu) = Sz'nri 

(u a upu~ ( ) = SzUanpn-f H (S a/ 3n 7 + S^np + Sp^n a ). (40) 

o 

E.g. the distribution functions (1161 and d!7i are uniaxial symmetric with respect to the direction of the magnetic field, n = H. 
Expressions for the equilibrium order parameters S? q are given in Eqs. d!9l and d20i for i = 1,2. 

Inserting Eqs. J40t into ( I35> and ( I36> and using Eq. ( 12 8 1 the hydrodynamic stress tensor is of the form assumed in the 
Ericksen-Leslie theory of nematic liquid crystals 01, 

T = a\(D : nn)nn + a2nN + a^Nn + a^D + a§nn ■ D + a$D ■ nn (41) 

where N = n — fixn is the corotational derivative of the director n and the Leslie coefficients a, are given by 

o c-2 o c2 

ai = 0, a 2 = -3t? s 0— -\- + l^aX^ 2 S\, a 3 = 3r/ s 0— -L- + 7^*7 A0 2 S 2 (42) 



a 4 = 2 Vo + 4(1 - d/^aX^Sf, a 5 = (ci - 13/2)2^ A0 2 5 2 , a 6 = (c x + l/2)2r ?s aA0 2 5 2 . (43) 
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A relation between die Leslie coefficients can be derived from Onsager's reciprocity relation, a 2 + a 3 = ag — as, which is 
known as Parodi's relation B Usll . Parodi's relation is readily verified from Eqs. (1421 and ( 143 \ . Note, that in the limit A — > 0, the 
result of Ref. I20I1 is recovered from Eqs. J4-1I . J42i and d43i . 

The balance equation for the director n can be derived from the moment equation (12 8 1 with the help of Eqs. J40i . 



with the molecular field H 



mol 



(1 — nn) ■ 
nksTSih and 71 



11 



mol 



— 71 AT — 72 D ■ n 



= 0, 



(44) 



a 3 - a 2 and 72 = a 6 ~ a 5 . 



A. Effective Field Approximation 

In the so-called Effective Field Approximation (EFA) @], a special family fe e of uniaxial distribution functions is considered 
that is obtained by replacing the magnetic field h in Eq. (I17t with an effective field £ e . Motivated by the good accuracy of the 
EFA for the non-interacting model @], we consider the EFA also for the present model. 

In more general terms, the EFA can be interpreted as the Quasi-Equilibrium Approximation where only the magnetization 
is kept as macroscopic variable (see e.g. fl5lo . Extremizing the free energy functional dlOi subject to the constraint of fixed 
normalization and fixed value of the first moment yields the quasi-equilibrium distribution fe e . Also in the present case, 
is obtained from the equilibrium distribution / eq if the magnetic field h in Eq. d!6> is replaced by an effective field £ e . Since 
the equilibrium distribution / cq is uniaxially symmetric around the magnetic field direction H, the distribution function fe e is 
uniaxially symmetric with respect to the direction of the effective field n = £ e , where £ e = £ e £ e and £ e denotes the norm of £ e . 
Consequently, the result of the previous section apply to the EFA. In particular, the moments (it) and (uu) within the EFA are 
given by Eq. J40i . where the scalar orientational order parameters are obtained from their equilibrium values by Sj = S* q (£ e ). 

1. Magnetization Equation 

From the moment equation (12 8 1 a closed equation for the magnetization M = M sat (it) can be derived within the EFA which 
reads 

M fix M = -—A - —A ■ MM + X 2 D -M + \ 3 D: MMM, (45) 

v\ v 2 

where A = kBT(£ e — h)/^i is the (dimensional) deviation of the effective field from the magnetic field. The coefficients i/j and 
Xi are defined as 

- = 3D rXL A^e), ± = -J*-B(£ e ), (46) 

V\ V 2 OT] s (p 



A 2 = -a x L A(Ce), A 3 = ^B(^). (47) 



The functions A(£ e ) and B(£ e ) are given by 



A(« - 2 -±^- = 1 - ^ + X.^iL^f L^ e) ) + lc 2 ^,m e ), (48) 
Wt \ - 52 fej j- (i + T <? U1 5 + ^(6) A , ^ 2 L 1 ^ e ).m e )-2L 2 {^)G' 2 ^ e ) 



and 



The magnetization equation (145 \ is a special case of Eq. (15) of Ref. 12111 which has been derived within a thermodynamic 
framework. The coefficients appearing in the magnetization equation, however, cannot be determined within the thermodynamic 
approach. For the special case A; = 0, Eq. ( 145 \ has been derived in Ref. [9] in linear order in A and <f>, within the EFA, starting 
from an A^-particle Fokker-Planck equation. The expression for A(£ e ) given by Eq. (I48i is identical to the result of Ref. 0] to 
first order in A, while the result of Ref. 01 for the coefficient B(£ e ), coincides with Eq. ( I49> only for A = 0. However, correcting 



Eq. (15) of Ref. ]9|] for the missing factor Li(£ e )/£ e B22II . also the results for the coefficient B(£ e ) agree with Eq. d49t . In 
Ref. IU5I1 . the magnetization equation ( 145 \ has been derived within a kinetic model of non-interacting, ferromagnetic colloidal 
particles with an ellipsoidal shape. Comparing Eq. 1471 to Eq. (77) of Ref. 11511, we notice that the transport coefficients Aj are of 
a similar form in both cases. For weak fields, the coefficient A 2 approaches a constant value A 2 (0) = — |crox L which is identical 
to the result of a non-interacting system of particles with shape factor B = — 4^oox L , B = (r 2 — 1)/ (r 2 + 1) where r is the 
axis ratio. Note, that B < corresponds to oblate, B > to prolate ellipsoidal particles. 
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2. Relaxation Times 

Analytical results for the magnetization dynamics ( 145 \ can be obtained for small deviations from the equilibrium values. To 
linear order in A and in the absence of velocity gradients, Eq. ( 145 \ becomes 

M=-4qA-4qA-M cc >M^, (50) 

where v cq denote the equilibrium values of the coefficients vi. Decomposing the off-equilibrium magnetization into components 
parallel and perpendicular to the magnetic field direction, M = Af" + 7Vf x , one finds from Eq. &50\ 

M = — ^r-M 1 - — \r(M^ - M cq ), (51) 
where the field-dependent relaxation times are defined by 

. = 3^(a) „ = 35T(a) 

D r a{2 + 5 2 cq (a)) ' 2D r (l - S° q {a)) ' V ; 

Linearization in the volume fraction (f> leads to 

r X =T^(1 + X ji + X 2 4>4) (53) 



where 



L\{a) || _ aL[(a) 

D r (a - Li(a)) ' T ° ~~ 2D r L 1 (a) 



(55) 



are the corresponding relaxation times in the non-interacting system and the functions tj- and t\ are defined in appendix|B] For 
the case of vanishing magnetic field, a — > 0, the above expressions coincide, t ± (0) = t"(0) = r rot (l + Xl/3). Note, that no 
contribution from C(A 2 ) remains in this limit. Thus, dipolar interactions lead to an increase of the zero-field relaxation time 
compared to the dilute suspension. In Fig. [3] we plot the relaxation times t 1 - and r" as a function of the magnetic field a for 
4> = 0.15 and A = 0, 1, 1.5, respectively. From Fig.|3]we notice that the transverse relaxation time is enhanced due to dipolar 
interactions for arbitrary values of the magnetic field. As the magnetic field increases, however, differences between r and 

decrease. The relaxation time parallel to the magnetic field is increased compared to rj only for small magnetic fields, while it 
is decreased in case of strong magnetic fields. For comparison, we included in Fig.|3]also the corresponding results of Ref. 01 . 
Note, however, that due to the dependence of t 1 - and Tq~ on the coefficient B, the results of Ref. 01 for the relaxation times are 
incorrect ll22ll . 

V. RESULTS FOR PLANE COUETTE FLOW 

To further illustrate the result (I41H43I . we consider in the sequel the special case of plane Couette flow, v(r) = (jy, 0, 0), 
where 7 denotes the constant shear rate. If the magnetization is oriented in the flow, gradient and vorticity direction, the 
Miesowicz shear viscosity rji = (0:3 + + a§)/2, r\2 = (— «2 + «4 + Oin)/2 and 773 = 0:4/2 is measured, respectively. From 
Eqs. (E2J and @3} we find 



Vi = Va + V? 



oc2 



2 + 5 : 



(56) 



where 77^° = |?7 S 0, is the Kronecker symbol, di = (ci/6 + 3)/6, d 2 = ( c i — 4)/6, and d 3 = (1 — Ci/3)/6. FromEq. ( I56l > 
we find -qi > i] 2 - The same inequality is found for a dilute suspension of ferromagnetic oblate ellipsoidal particles Hal , while 
771 < 1J2 is valid for oblate, chain-like particles. 

If the magnetic field is sufficiently strong to fully orient the magnetic moments, Si — > 1, the Miesowicz shear viscosities 
approach their asymptotic values rji — > 77°°. From Eq. d56l one finds 



Vt° = Vo + t?~(1 - « <l3 + di^xj. (57) 
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Since 77™ is the maximum viscosity increase for A = 0, we find from Eq. i57l that dipolar interactions increase the value of 
77^° while 77!° is decreased (for C\ < 24) compared to 77^°. Eq. J57l > further predicts the inequality rj^° > 77^° > 77^°, which is 
in qualitative agreement with the results of nonequilibrium molecular dynamics simulations on a magnetically saturated model- 
ferrofiuid in the weakly interaction regime A < 1 l23ll . 

For small shear rates, the order parameters Si in Eq. d56t may be replaced by their equilibrium values. As has been pointed 
out in Ref. fl5ll . the assumption of uniaxial symmetry reduces to the EFA in this limit. Inserting Eqs. d!9l > and J20t into Eq. d56l > 
and expanding in the small quantities A and <f> gives 



771,2 (a) = rjo + r)7 



ah\(a) 
a — L\(a) 



1 + xAti(a) + L\ (a) + di,a(l - ^-)a} + A0 2 {^(a) + ^M} 



(58) 



and 773(a) = 770 + d^ax^L^a) . The functions tf~ are defined in the appendixlBl For a — > 00, the result J57t is recovered 
from Eq. ( I58> . In the opposite limit of weak magnetic fields, a — > 0, the viscosities 77, increase quadratically with a, 



rn(a) = 770 + 7 T7s0 1 - &,3 + -* L (1 - + ha) a 1 + 0(a A ) 



(59) 



We now consider the special case where the magnetic field and the magnetization are oriented in the plane of shear. In this 
case, we introduce the angles of the flow direction with the magnetic field 7? and with the magnetization 9, respectively, such that 
H = (cos 7?, sin 7?, 0) and n = (cos 9, sin 8, 0). Taking the vector product of the moment equation (12 St with n, the alignment 
angle 9 of the magnetization in the stationary state can be found, 



sm(9 - 1?) + Mn 



35i 



aox^Si cos(26») 



0. 



(60) 



In Eq. J60l >. we have introduced the Mason number Mn = T ro t7/o;, which measures the relative strength of the flow compared 
to the magnetic field. From Eq. J60t we notice that a perfect alignment of the magnetization with the magnetic field, 9 = -d, 
occurs for Mn = 0. For Mn > 0, the direction of the magnetization does not coincide with the magnetic field direction due to 
the effect of the flow field. With the ansatz 6 = J2k Mnfe #fc, where 9 = 7? denotes the value of 9 for Mn = 0, Eq. (|60j can be 
solved recursively for 9. The first order result reads 



occq 

—l^ + ^S^cosW) 



(61) 



Due to dipolar interactions, the present model also predicts normal stress differences. The first normal stress coefficient N\ 
for the plane shear flow v = (jy, 0, 0) is defined as Ni = —(T xx — T yy ), while the second normal stress coefficient is defined 
as N2 = —(Tyy — T zz ). From Eqs. (141 1 — d43i . we find that iVj = 14ipiri s a\(t) 2 Sfjn x n y if the magnetic field is oriented in the 
plane of shear, where ijji = 1, 1P2 = — (cj + l/2)/2 and n x , n y denote Cartesian components of the director n. Note, that 1P2 is 
the ratio N2/N1. For weak velocity gradients, the normal stress coefficients take the form 



Ni = 



Uip i r) a a\4> 2 Sf'y [sin 7? cos 7? + Mn0i cos(2tJ»)] . 



(62) 



From Eq. d62i we observe that Ni cx 7 for weak velocity gradients if the magnetic field is oriented either in flow or in gradient 



direction, otherwise we find Ni cx 7. In any case, the normal stress differences vary as Ni cx a 
weak, r rot 7 < 1. 

In a pipe flow, the so-called McTague Jll y[ viscosity coefficients 7711 = 7/1 and 77^ = (772 
magnetization is oriented in flow and perpendicular to the flow direction, respectively. From Eq. 

A^ 



V\\ - Vo 



2 1 + 2d!a Xl 



S 2 



for a -c 1 if the flow field is 

+ 773 )/2 are measured, if the 
l56l we find to lowest order in 



(63) 



From Eq. d63l > we observe that dipolar interactions lead to a ratio of the McTague viscosities greater than two, the ratio being an 
increasing function of the magnetic field. It has been emphasized several times that McTague's experimental results 0] can be 
fitted quite nicely by the predictions of kinetic model introduced in Ref. p4:]. However, the ferrofluid used in these experiments 
cannot be considered dilute and the parameters of the fit (the volume fraction <j> and the magnetic moment of the colloidal 
particles /j) do not agree with independent estimates 0]. In addition, the ratio d63i found in the experiments 1 1] is not constant 
and becomes approximately 2.3 for strong magnetic fields. This result cannot be explained within the model |4] which predicts 
a value of 2 independent of the magnetic field. In view of the present results, a value greater than two for the ratio Eq. (1631 
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is expected due to magnetic dipole-dipole interactions. It would be very interesting, to compare the prediction of Eq. d63i to 
experimental results on a semi-dilute ferrofluid with A < 1 where the present model applies. 

The maximum relative viscosity increase if the magnetization is oriented in flow direction, Ajj° = (r/j — r]o)/r/o, is given by 

+ O(0 3 ), (64) 

where Ajj° = 30/2 corresponds to the classical result of Ref. 0|. Note, that the factor 5/2 in Eq. J64> arises from the expansion 
of ?7o in terms of <f>. Experimental measurements of Ay have been performed in Ref. @] with a capillary viscosimeter for a 
magnetite based ferrofluid with volume fractions <j> < 0.2. The interaction parameter A, Eq. (J3jl, was estimated in Ref. |7] to 
be A = 0.2, justifying the application of the present model to their system. Fig. |4]shows Ajj° (cf>) according to the experimental 
results of Ref. 0]. From Fig.|4]it is seen that Ajj° raises stronger than linearly with </> for <j> > 0.1 and that Eq. J64i describes 
the experimental data accurately up to volume fractions <f> < 0.15. A value of d\o — 2.75 was used in Fig. [4] For a soft sphere 
potential we obtained c\ = 12 in Sec. [HI] so that d\o = 2.75 corresponds to a as 2.3 which is within the expected range of 
values. For higher concentrations, the present approach needs to be extended by including higher order terms in the expansion 
©. 

Figure |5] shows the dimensionless viscosity increase A|| = (j]i — r]o)/r]o as a function of a for various values of 0. The 
experimental results for A|| of Ref. 0] are shown together with Eq. fl58i for i = 1. As in Fig.0] a value dicr = 2.75 has been 
chosen. From Fig. [5] we notice that Eq. (1581 is able to describe the experimental data well for volume fractions <f> < 0.1 while 
strong deviations appear for <fi = 0.17. It has been mentioned in Ref. |7] that their experimental results can be fitted to the 
non-interacting model, A = 0, with the help of an effective volume fraction <f> c g > <j>. Within the present model, the increase of 
the effective volume fraction is explained by the magnetic dipole-dipole interaction. 

VI. CONCLUSION 

In the present work, we have proposed a kinetic model of dilute, weakly interacting ferrofluids that extends the classical kinetic 
model of dilute ferrofluids by the incorporation of dipolar interactions. Our model predicts several extensions compared 
to the classical kinetic model which are in qualitative agreement with experimental results, such as the presence of normal 
stress differences, enhanced magnetoviscous effect, modified anisotropy of viscosity and the dependence of viscosity on the 
hydrodynamic volume fraction and the symmetric velocity gradient. For a quantitative comparison with experimental results, two 
additional parameters have to be specified which are absent in the non-interacting model 0: the dipolar interaction parameter 
A and the ratio of the rotational over the translational relaxation time of the colloidal particles a. While A is defined by Eq. Q 
and also tabulated for several ferrofluids, determining a is not straightforward. As has been discussed above, estimations of a 
can give information only about the order of magnitude of a, rather than a certain value. On the other hand, the value of a can 
be inferred from measurements, e.g. of the maximum viscosity increase from Eq. (I57> . Following this route, the experimental 
results of Ref. are described quantitatively by the present model with a reasonable value of a. 
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APPENDIX A 

The functions Gi(a s ) are defined by Gi(a s ) = Gi[f cq ], where the functionals G,[/] for i = 2, 3, 4 are given by Eqs. (II 1> - 
(1131 and / eq is defined in Eq. d!6i . Evaluating the functionals Gi[f] with the equilibrium distribution function / oq one obtains 

G 2 (x) = -i (L 2 (x) 2 + 5) (Al) 



G 3 (x) = --!-(L 3 (x) 2 -21L 1 (x) 2 ) 



(A2) 
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Ga{x) = -i- (L 4 (x) 2 + 20L 2 (x) 2 + 49) . (A3) 
3675 v ' 



The functions G,; (x), defined in Eqs. JA1HA3> . are monotonously increasing functions of x. These functions depend only weakly 
on x and have the following expansion for x — > and x — > oo: 



t + 3375 x 4 + 0(x 5 ) for x -> 
1(1 - a:" 1 ) + e>(x~ 2 ) forx^oo 



G 2 (X) = ^ 8 V 3375- ^ (A4) 



16 „4 



x 4 + 0(x 5 ) for x -> 



G A {x) = { 8 VV 5 -nln - w (A6) 
APPENDIX B 

The contribution of dipolar interactions to the transverse and parallel relaxation times are described by the functions 

ti(x) = L[(x) - Ll{ T X ] J L^xf - L 2 (x)) (Bl) 

h[X) ~L 1 (x) 3(x-L x (x)) { } 



t „ Mx) + 2Mx) 
x L x (x) 



2 r / m \ \ L 2(x) 



L 1 (x)(L 1 (x)"- J L 2 (x))- 



(B3) 



II, . G" 2 '(x) xJ'(x) 

The functions tj- are positive, while functions t| are not sign-definite. In the limit x — -> oo we find tf- -> and tf -> 0, while 
for x — > 0, the following asymptotic behavior is obtained: tj-(x) = 1/3 - 4x 2 /45, t 2 {x) = Ax 2 /375, t\ (x) = 1/3 - 7x 2 /45, 
and if (x) = 56x 2 /1125. 
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Langevin parameter 

FIG. 1: Equilibrium magnetization as a function of the Langevin parameter for volume fraction 4> = 0.157 and A = 1. Symbols are the result 
of molecular dynamics simulations presented in Ref. fI3l . Dashed line is the Langevin function L\(a), while the solid line corresponds to 
Eq. 1 1 9i . where the infinite sum has been truncated at k = 4. The dotted line is the result of the approximation L\ (a + x L Li (a)). 




FIG. 2: Shear viscosity tjo/^s of the isotropic suspension (H — 0) as a function of the volume fraction <f). Dashed an full lines are the 
theoretical predictions rjo/^s = 1 + 2.5(f) + b(f> 2 with 6 = (Einstein) and b = 6.2 (Batchelor), respectively. 
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FIG. 3: a) Reduced transverse, r ± /r ro t , and b) parallel, r"/r ro t , relaxation times as a function of the applied magnetic field a. The volume 
fraction was chosen as cj> — 0.15. Curve 1 corresponds to A — 0, curve 2 to A = 1, and curve 3 to A = 1.5. Solid lines are the result of 
Eqs. J53I and 1541 . dashed lines are the result of Ref. 0]. 
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volume fraction 



FIG. 4: The maximum relative viscosity increase Ajj° as a function of the hydrodynamic volume fraction cj>- Circles represent the experimental 
results from Ref. 01, solid line the prediction of the present model, Eq. 1641 . for dia = 2.75 and dashed line is the result of Ref. 01. 
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FIG. 5: Relative viscosity increase Ay as a function of the Langevin parameter. Circles, squares and diamonds represent the experimental 
results of Ref. 01 for volume fractions <j> — 0.17, 0.11 and 0.06, respectively. Solid lines are the corresponding predictions of the present 
model, Eq. {58} with dia = 2.75. 



